Figure 4

Immunoglobulin plots

Heatmap of significant genes per cluster

Load data and libraries

##################
# LOAD LIBRARIES #
##################
library(tidyverse)
library(Seurat)
library(SeuratObject)
library(tidyseurat)
library(cowplot)
library(patchwork)
library(openxlsx)

source("../../bin/spatial_visualization.R")
source("../../bin/plotting_functions.R")

#########
# PATHS #
#########
input_dir <- "../../results/06_DGE_condition_st_data/"
result_dir <- "./Figures/05/"
if( isFALSE(dir.exists(result_dir)) ) { dir.create(result_dir,recursive = TRUE) }
epi_clus <- "^5$|^6$|^7|^8" # res 0.7

ord <- c("Superficial", "Upper IM", "Lower IM", "Basal","1","4","0","3","2","9","10","11","12")
ord1 <- c("5", "6", "7", "8","1","4","0","3","2","9","10","11","12")
sample_id <- c("P020", "P045", "P050", "P057",
               "P008", "P031", "P044","P080", "P026", "P105", 
               "P001", "P004", "P014", "P018", "P087", "P118",
               "P021", "P024", "P067", "P081", "P117" ) %>% set_names()

#############
# LOAD DATA #
#############
DATA <- readRDS(paste0("../../results/03_clustering_st_data/","seuratObj_clustered.RDS"))
DEGs_table <- read_csv(paste0(input_dir,"DGEs_condition_wilcox.0.7.csv"))
pseudo_DEGs <- readRDS(paste0(input_dir,"Pseudobulk_across_DEGs.RDS"))

D <- DEGs_table %>%
  filter(., p_val_adj < 0.05) %>%
  filter(., !(grepl("0|^3|9|10|11|12", .$subgroup))) %>%
  #filter(., cluster == "L4") %>%
  mutate(sp_annot = ifelse(grepl("\\d", .$layers), "SubMuc", "epi")) %>%
  split(~sp_annot) 

# NOT USED
#################
# PLOT FUNCTION #
#################
morf_DEGs_plot <- function(genes = genes_epi, clus="^5$|^6$|^7|^8"){
  df <- DATA %>%
    mutate(., FetchData(., vars = c(names(genes)), slot = "counts")) %>%
    filter(., grepl(clus, .$Clusters)) %>%
    as_tibble() %>%
    select(1:5, layers, all_of(names(genes))) %>%
    pivot_longer(cols = any_of(names(genes)), names_to = "gene", values_to = "values") %>%
    # filter(gene == "LGALS7") %>%
    group_by( groups, gene, layers) %>%
    summarize(sum_counts = sum(values), .groups="drop") %>%
    mutate("Sum counts(log10)" = log10(.$sum_counts)) %>%
    arrange(`Sum counts(log10)`) %>%
    mutate(type = genes[.$gene])
  
  absmax <- function(x) { x[which.max( abs(x) )]}
  ord <<- df %>%
    #pivot_wider(id_cols = -`Sum counts(log10)`,names_from = "groups", values_from = "sum_counts") %>%
    pivot_wider(id_cols = -sum_counts, names_from = "groups", values_from = `Sum counts(log10)`) %>%
    mutate(diff_L1_L4 = L4-L1, 
           diff_L2_L4 = L4-L2, 
           diff_L3_L4 = L4-L3, .by = "gene") %>%
    rowwise() %>% 
    # Identifies the larges value of eah row among the L1-L4 columns:
    mutate(Regulation = sum(c_across(starts_with("diff")) ),
           max = paste0(names(.[4:7])[c_across(starts_with("L",ignore.case=F)) == 
                                        absmax(c_across(starts_with("L",ignore.case=F)))], collapse = '_') ) %>%
    ungroup() %>%
    mutate(Regulation = ifelse(.$Regulation < 0, "DOWN", "UP"),
           max = str_extract(.$max, "L\\d"))
  
  m <- ifelse(clus == "^5$|^6$|^7|^8", "Epithelial", "Submucosal")
  #write.xlsx(ord, paste0(result_dir, "Selected_DEGs_FIG2-3_", m, ".xlsx"))
  #####################
  # TOP DEGs PLOTTING #
  #####################
  g_ord <- arrange(ord, diff_L1_L4) %>% pull("gene") %>% unique()
  
  col <- c("L1"="#FF7F00", "L2"="#FED9A6", "L3"="#6A51A3", "L4"="#9E9AC8" ) # "#FFFFB3", "#BEBADA","#8DD3C7", "#FB8072",
  col <- c("L1"="#56B4E9","L2"="#009E73","L3"="#CC79A7", "L4"="#FC8D62")
  
  (B <- left_join(df, select(ord, Regulation, max, gene), by="gene") %>%
      #mutate(., gene = factor(gene, levels=g_ord)) %>%
      ggplot2::ggplot(data=., aes(x=fct_reorder2(gene, Regulation,`Sum counts(log10)`), y=`Sum counts(log10)`)) +
    geom_point(aes(col=groups), size = 2) +
    scale_colour_manual(values = col) +
    facet_grid(cols=vars(), rows = vars(layers), scales = "free_x", space = "free_x") +
    #facet_wrap("max", ncol = 2, strip.position = "top", scales = "free_x", shrink = F) +
    theme_minimal() +
    theme(legend.position = "bottom",
          plot.margin = unit(c(5,-1,5,1),units = "pt"),
          legend.margin = margin(-15,2,-8,2), # moves the legend box
          legend.title = element_blank(),
          legend.text = element_text(size = 8),
          legend.spacing.x = unit(4, 'pt'),
          axis.title = element_text(size=8),
          axis.text.x = element_text(size=8, angle=45, hjust=1, color="black"),
          axis.title.x = element_blank(),
          #text = element_text(size = 10),
    ) )
}
geneid  <- c("IGHG4", "IGLC2", "IGHG1", "IGHA2", "IGHA1")
genes <- DEGs_table %>% filter(gene %in% geneid) %>% filter(grepl("^1$|Basal", .$layers)) %>% distinct(gene, layers) %>% deframe()
(Immuno <- morf_DEGs_plot(genes = genes, clus="^1$|8"))

###############
# VIOLIN PLOT #
###############
# https://stackoverflow.com/questions/35717353/split-violin-plot-with-ggplot2

feature <- c("IGHG4", "IGLC2", "IGHG1", "IGHA2", "IGHA1", "CD81", "IGKC", "JCHAIN")
feature <- c("IGHA1","IGHG1", "CD81", "IGKC", "JCHAIN")
col <- c("L1"="#56B4E9","L2"="#009E73","L3"="#CC79A7", "L4"="#FC8D62")
facet <- "layers"
group.by <- "groups"
ncol <- 1

DEGs_df <- DATA %>%
  #filter(!(grepl("\\d", .$layers))) %>%
  filter((grepl("Basal|^1$", .$layers))) %>%
  filter((grepl("L1|L2", .$groups))) %>%
  mutate(., FetchData(., vars = c(feature)) ) %>%
  #mutate(across(any_of(feature), na_if, 0)) %>% # NB! removes zero values 
  select(any_of(c(feature, facet, group.by, "layers"))) %>%
  pivot_longer(cols = any_of(feature), names_to = "feature", values_to = "val") %>%
  mutate(comb = paste0(.$groups, "_",.$layers)) %>%
  mutate(groups = factor(.$groups, levels = sort(unique(.$groups))))

# groups on y-axis, genes on x axis
# not used:
violin_mean_plot <- function(DEGs_df){
  m <- max(na.omit(DEGs_df[["val"]]))/1 # try e.g 2
  DEGs_df %>%
    ggplot(aes(y=.data[[group.by]], x=.data[["val"]], col=.data[[group.by]])) +
    geom_violin() + # {if(length(feature == 1)) ggtitle(feature)} +
    #geom_boxplot(fill = "transparent", notch = T, outlier.size = .5, outlier.colour = "red") +
    geom_jitter(width = 0.2, height = 0.3, alpha = 0.2, size=.1) +
    stat_summary(aes(x = val), fun=median, geom="crossbar", linewidth=.3, color="black") +
    scale_color_manual(values = col) +
    my_theme + NoLegend() + #xlim(c(0, m)) +
    theme(text = element_text(size = 10)) +
    #       axis.text.x = element_text(angle = 30, hjust=1),
    #       plot.title = element_text(hjust = 0.5),
    #       axis.title.y = element_blank(),
    #       axis.title.x = element_blank()) +
    #facet_wrap(~.data[[facet]], ncol = ncol, strip.position = c("right", "bottom", "left", "top"))
    facet_grid(rows = vars(layers), cols = vars(feature)) #, strip.position = c("right", "bottom", "left", "top"))
}

# dev.new(width=6.7, height=3, noRStudioGD = TRUE) 
A <- violin_mean_plot(DEGs_df)
# ggsave("./Figures/04/Fig_viol.png", A, width = 6.7, height = 3) #, dpi = 300

# groups on x-axis, genes on x axis
# used for figure 4
violin_mean_plot <- function(DEGs_df){
  m <- max(na.omit(DEGs_df[["val"]]))/1 # try e.g 2
  DEGs_df %>%
    ggplot(aes(x=.data[[group.by]], y=.data[["val"]], col=.data[[group.by]])) +
    geom_violin() + # {if(length(feature == 1)) ggtitle(feature)} +
    stat_summary(aes(y = val), fun=median, geom="crossbar", linewidth=.3, color="black") +
    #geom_boxplot(fill = "transparent", notch = T, outlier.size = .5, outlier.colour = "red") +
    geom_jitter(width = 0.3, height = 0.2, alpha = 0.2, size=.01) +
    scale_color_manual(values = col) + ylab("Expression") +
    my_theme + NoLegend() + #xlim(c(0, m)) +
    theme(text = element_text(size = 12),
          axis.line = element_blank(),
          panel.border = element_rect(colour = "gray") ) +
    #       axis.text.x = element_text(angle = 30, hjust=1),
    #       plot.title = element_text(hjust = 0.5),
    #       axis.title.y = element_blank(),
    #       axis.title.x = element_blank()) +
    #facet_wrap(~.data[[facet]], ncol = ncol, strip.position = c("right", "bottom", "left", "top"))
    facet_grid(rows = vars(layers), cols = vars(feature)) #, strip.position = c("right", "bottom", "left", "top"))
}

# dev.new(width=3.3, height=3, noRStudioGD = TRUE) 
A <- violin_mean_plot(DEGs_df)
# ggsave("./Figures/04/Fig_L1_L2_viol.png", A, width = 4, height = 4, dpi = 500) #, dpi = 300

######################
# PLOTTING ON TISSUE #
######################
col = c("#EFEDF5", "#DADAEB", "#BCBDDC", "#9E9AC8", "#807DBA", "#6A51A3", "#54278F", "#3F007D")
geneid <- c("IGHG4", "IGLC2", "IGHG1", "IGHA2", "IGHA1", "CD81", "IGKC", "JCHAIN")

# ggsave("./Figures/02&03/Tissue_Immuno.png", p, width = 12.5, height = 12.5, bg = "white") #, dpi = 300
# dev.new(height=3, width=7, noRStudioGD = TRUE)
plots_all <- geneid %>%
  map(., ~plot_spatial.fun(
      DATA,
      sampleid = sample_id,
      save_space = T,
      colors = col, 
      geneid = .x,
      zoom = "zoom",
      ncol = 4,
      img_alpha = 0,
      point_size = .5)
    )

# dev.new(width=8, height=8, noRStudioGD = TRUE)
plots_all[[1]]

################
# SAVE RESULTS #
################
pdf(file=paste0("./Figures/04/Tissue_Immuno_all.pdf"), 
    #width = 8.5, height = 10 # Epi
    width = 8, height = 8 # Sub
)
plots_all

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

dev.off()

## quartz_off_screen 
##                 2

######################
# PLOTTING ON TISSUE #
######################
sample_id_ <- c("P057","P044")


plots_rep <- geneid %>%
  map(., ~plot_spatial.fun(
      DATA,
      sampleid = sample_id_,
      save_space = F,
      colors = col, 
      geneid = .x,
      zoom = "zoom",
      ncol = 4,
      img_alpha = 0,
      point_size = .5)
    )
# dev.new(width=3.3, height=2, noRStudioGD = TRUE)
plots_rep[[1]]

B <- plot_grid(plots_rep[[8]], plots_rep[[5]], ncol = 1)
# ggsave("./Figures/04/on_tissue.png", B, width = 4, height = 3.3, bg = "white", dpi = 500) #, dpi = 300

#########################
# COMBINE PANEL A AND B #
#########################
# dev.new(width=6.7, height=3, noRStudioGD = TRUE) 
(FIGURE4 <- plot_grid(A, B, rel_widths = c(1, 1))) # labels = c("a", "b")

# ggsave("./Figures/04/Figure_4.png", p, width = 6.7, height = 3.3, bg = "white") #, dpi = 300